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The  next  generation  of  electronic  devices  will  be  developed  at  the  nanoscale  and  molecular  level, 
where  quantum  mechanical  effects  are  observed.  These  effects  must  be  accounted  for  in  the  design 
process  for  such  small  devices.  One  prototypical  nanoscale  semiconductor  device  under 
investigation  is  a  resonant  tunneling  diode  (RTD).  Scientists  are  hopeful  the  quantum  tunneling 
effects  present  in  an  RTD  can  be  exploited  to  induce  and  sustain  THz  frequency  current  oscillations. 
To  simulate  the  electron  transport  within  the  RTD,  the  Wigner-Poisson  equations  are  used.  These 
equations  describe  the  time  evolution  of  the  electrons’  distribution  within  the  device.  In  this  paper, 
this  model  and  a  parameter  study  using  this  model  will  be  presented.  The  parameter  study  involves 
calculating  the  steady-state  current  output  from  the  RTD  as  a  function  of  an  applied  voltage  drop 
across  the  RTD  and  also  calculating  the  stability  of  that  solution.  To  implement  the  parameter  study, 
the  computational  model  was  connected  to  LOCA  (Library  of  Continuation  Algorithms),  a  part  of 
Sandia  National  Laboratories  parallel  solver  project,  Trilinos.  Numerical  results  will  be  presented. 
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1.  Introduction 

In  the  next  few  decades,  if  the  trend  predicted  by  Moore’s  Law1  is  to  continue, 
semiconductor  devices  will  need  to  be  scaled  down  to  the  nanoscale.  The  principal 
hindrance  in  the  design  of  these  small  devices  will  be  that  quantum  physics,  instead  of 
classical  physics,  will  dictate  their  operation.  Physicists  and  engineers  do  not  yet  fully 
understand  what  new  features  will  be  present  due  to  quantum  effects.  This  lack  of 
knowledge  motivates  deriving  an  accurate  model  of  electron  transport  at  the  quantum 
level.  One  such  model,  the  Wigner-Poisson  equations2,  is  utilized  in  this  paper  for 
simulating  electron  transport  in  a  prototypical  nanoscale  semiconductor  device,  the 
resonant  tunneling  diode  (RTD). 

The  RTD  has  been  extensively  researched  in  the  past  two  decades3' 4’  \  In  the  1980s, 
these  devices  were  placed  in  circuit  configurations  as  gain  elements  to  induce  a  high- 
frequency  current  oscillation6.  While  high-frequency  oscillations  developed  (high- 
frequency  here  means  THz),  oscillations  at  lower  frequency  modes  also  developed,  and 
too  much  power  was  lost  at  these  lower  frequency  oscillations.  To  prohibit  these  lower 
frequency  oscillations  from  developing,  the  device  size  had  to  be  further  scaled  down,  but 
this  also  led  to  reduced  power  output  from  the  high  frequency  oscillation7.  Ongoing 
research  with  the  RTD  involves  removing  it  from  the  circuit  to  see  if  some  intrinsic 
mechanism  can  be  exploited  to  created  THz  frequency  oscillations. 

Figure  1  presents  a  schematic  of  an  RTD  and  the  diagram  of  the  electric  potential 
within  the  RTD  when  a  large  voltage  drop  V  is  applied  to  the  RTD.  An  RTD  is  made  by 
growing  alternating  layers  of  semiconductor  materials,  denoted  by  material  I  and  material 
II  in  Figure  1.  The  difference  between  the  materials  is  that  material  11  has  a  larger  band- 
gap  than  material  I.  This  difference  causes  the  discontinuous  jumps  in  the  electric 
potential  present  in  Figure  1.  For  our  simulations,  material  I  is  gallium  arsenide  (GaAs), 
and  material  II  is  aluminum  gallium  arsenide  (AlGaAs). 

At  the  left  and  right  ends  of  the  RTD,  there  is  doping,  where  atoms  that  have  more 
electrons  than  the  semiconductor  is  placed  in  the  device  to  increase  the  number  of 
electrons  present  in  the  device.  The  voltage  drop  across  the  device  drives  electrons  to  the 
right  side  of  the  RTD,  creating  a  current.  Since  in  quantum  physics,  electrons  are  treated 
as  waves  instead  of  particles,  when  they  encounter  the  potential  barriers,  they  have  the 
ability  to  tunnel  through  them  and  reach  the  right  side  of  the  RTD.  This  quantum 
tunneling  effect  is  the  basis  of  the  RTD’s  operation,  and  scientist  hope  this  mechanism 
can  be  understood  and  controlled  by  using  the  Wigner-Poisson  simulations  in  order  to 
create  and  sustain  THz  frequency  current  oscillation  within  the  RTD. 

2.  The  Wigner-Poisson  Equations 

The  Wigner-Poisson  equations  describe  the  time-evolution  of  the  electron  distribution 
within  the  RTD.  The  equations  consist  of  a  nonlinear  integro-partial  differential  equation 
(IPDE)  coupled  with  Poisson’s  equation  to  determine  the  electrostatic  potential  within  the 
device.  The  IPDE  is  given  by 
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Fig.  1 .  Diagram  of  an  RTD  and  the  electric  potential  within  the  RTD 


df/dt  =  W(f )  =  Kf  +  P(f)  +  S(f) .  (1) 

Here,/=/  ( x ,  k,  t)  is  the  electrons’  distribution  within  the  RTD.  It  is  a  function  of  the 
electron’s  position  x,  where  x  is  in  [0,  Z]  ( L  being  the  length  of  the  RTD),  electron’s 
momentum  k.  where  k  is  in  (-oo,  oo),  and  time  t>  0.  Given  an  initial  electron  distribution, 
Eq.  (1)  is  used  to  track  how  the  distribution  changes  in  time.  This  type  of  simulation  is 
discussed  in  the  time-dependent  simulation  section. 

The  second  type  of  simulation  we  perform  involves  directly  finding  steady-state 
electron  distributions.  These  distributions  will  satisfy  the  equation 

W(f)  =  Kf  +  P(f)  +  S(f)  =  0.  (2) 

This  type  of  simulation  is  discussed  in  the  steady-state  simulation  section. 

The  first  term  in  the  equation,  Kf,  is  a  linear  term  and  represents  the  kinetic  energy 
effects  on  the  distribution.  It  is  given  by 

Kf  =  -  hk/lmn  df  jdx .  (3) 

In  this  term,  h  is  Planck’s  constant,  and  m  is  the  electron’s  effective  mass. 

The  second  term  in  the  equation,  P  (/),  represents  the  potential  energy  effects  on  the 
distribution  and  is  the  nonlinear  term  in  the  equation.  It  is  given  by 
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/•OO 

P(f)  =  -  4 I/aJ  /(x,  k'  )T(x,  k  -  k'  )dk' .  (4) 

The  kernel  function,  T  (x,  k ),  which/is  multiplied  against  in  this  term  is 

T (x,  k )  =  f  [U (x  +  y)  —  U(x—  v)]  sin(2 yk)dy  .  (5) 

J  0 

Here,  U  (x)  is  the  electric  potential  within  the  device  (see  below),  and  Lc  is  the  correlation 
length.  The  term  is  nonlinear  in / since  U (x)  depends  on/ through  Poisson’s  equation. 

The  final  term  in  the  equation,  S  (/),  includes  scattering  effects  in  the  electron 
distribution,  and  its  formula  is 

s if)  =  1/ r[(J  x  f(x,  k)dk / j  f0  (x,  k)dk)f  (x,  k )  -f(x,k)\.  (6) 

In  this  term,  x  is  the  relaxation  time,  and  fo  ( x ,  k)  is  the  equilibrium  electron  distribution. 
This  is  the  solution  to  Eq.  (2)  when  there  is  no  voltage  drop  applied  to  the  device,  (i.e.  V 
=  0  in  Eq.  (10)  below). 

The  boundary  conditions  on/  are  imposed  at  the  incoming  boundaries.  For  the  left 
side  of  the  device  ( x  =  0),  this  means  imposing  boundary  conditions  for  k  >  0  since  for  k 
>0,  electrons  are  moving  from  left  to  right.  This  formula  is  given  by 

/(0,A:)  =  4  mnkBT  jh2  ln(l  +  exp[l/  k  BT(h2  k2  /8k2  m*  -//0)])-  (7) 

For  the  right  side  of  the  device  {x  =  L),  this  means  imposing  boundary  conditions  for  k  < 
0  since  for  k<  0,  electrons  are  moving  from  right  to  left.  This  formula  is  given  by 

f{L,k )  =  4 mnkBT jh2  ln(l  +  exp[l/£Br(7r£2/8;r2  rn  -  jUL)]) .  (8) 

In  the  boundary  condition  formulae,  kB  is  Boltzmann’s  constant,  T  is  the  temperature,  p0 
is  the  Fermi  energy  at  x  =  0,  and  |i/  is  the  Fermi  energy  at  x  =  L.  Figure  2  gives  a 
graphical  representation  of  the  boundary  conditions  in  (x,  k)  space. 

The  electrostatic  potential,  u  (x),  is  the  solution  to  Poisson’s  equation 

d2uldx2  =  q2  ls\NAx)  —  \l27t\  f(x,k)dk\.  (9) 

J-oo 

Here,  cj  is  the  charge  of  an  electron,  e  is  the  dielectric  permittivity,  and  Nj  (x)  is  the 
doping  profde.  The  boundary  conditions  on  u  (x)  are  where  the  voltage  drop  V  comes 
into  the  equations.  We  have 

u  (0)  =  0,  u  (L)  =  -V  (10) 

where  V>  0.  Once  u  (x)  is  solved,  then  the  electric  potential  is  calculated  by 

U  (x)  =  u  (x)  +  Ac  (x). 


(11) 


Simulating  Nanoscale  Semiconductor  Devices  5 


x=0 


(Adapted  from  Frensley,  Phys.  Rev.  B,  Vol.  36,  1987) 


x=L 


Doped  Doped 

Region  Region 


Fig.  2.  Boundary  conditions  on  the  electron  distribution  f  (x,  k ,  t) 


Here,  Ac  (x)  defines  the  discontinuous  jumps  in  the  electric  potential  created  by  the 
heterojunction  of  the  two  different  types  of  semiconductor  material. 

Finally,  to  calculate  the  current  density.  /  (x,  t ),  we  use  the  formula 

j(x,t)  =  h/lTim*  f  kf(x,k,t)dk.  (12) 

J- 00 


3.  Discretization 

To  numerically  solve  the  Wigner-Poisson  equations,  we  implement  a  finite  difference 
method.  This  discretization  converts  the  infinite-dimensional  IPDE  problem  into  a  finite- 
dimensional  nonlinear  ODE  problem  for  an  approximation  of f  at  the  grid  points  in  the  (x, 
k)  domain.  The  x-domain  of  [0,  L]  is  discretized  with  Nx  equally  grid  points  x,  =  (/'  - 

1  )Ax,  i=  1,2,...,  Nx  and  Ar  =  L/(NX  -  1).  The  ^-domain  of  (-oo,  oo)  is  truncated  to  (-Kmax, 
Kmax)  with  Nk  equally  spaced  grid  points  kj  =  (2 \j  -  Nk  -  l)Ak/2,j  =  1,2,  ...,  Nk  and  A k  = 

2  *Kmax/Nk. 

For  the  kinetic  energy  term,  Kf,  we  use  an  upwind  differencing  scheme  to 
approximate  the  spatial  derivative  term.  For  kj  >  0,  this  approximation  is  given  by 

Kf  *  -  hkj  / 2 mn  (-  3  ftj  +  4/MJ  -  ft_2J  /2Ax)  (13) 

and  for  kj  <  0,  this  approximation  is  given  by 
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Kf«-  hkt  jlmn  (if,  -  4  fMJ  +  ft2J /2to) 


(14) 


For  the  integrals  in  the  P  (/)  and  S  if)  terms  are  approximated  with  a  quadrature 
formula.  The  approximation  for  P  if)  is  given  by 


P(f)*-4/hYjfijT(xl,kJ-kr)Ak. 


(15) 


i= i 

The  approximation  to  the  kernel  function  T  (. x ,  k-k’)  is 


T(x„kj  -kr)  «  Yj[U(xi+xi,)-U(xi  - xv)\sm(2xv(kj  -kf))Ax.  (16) 


For  the  scattering  term,  the  quadrature  formula  gives 


W-mtfJtfMJMM.*,)-/,)-  on 

7=1  /  7=1 

Finally,  we  use  a  standard  three-point  center  differencing  scheme  to  numerically  solve 
Poisson’s  equation.  This  gives  for  i  =  2,  3,  . . .,  A^-1 

Nk 

m(x,_, ) - 2u{xi )  +  u(xM )/ Ax2  =q1/£[Nd(xi)-\/27rYJfiJ]  (18) 

7=1 


with  boundary  conditions  u(xi)  =  0  and  u(xNx)  =  -  V. 


4.  Time-Dependent  Simulation 

In  the  time-dependent  simulation,  the  time-evolution  of  the  electron  distribution  /  is 
tracked  from  the  ODE  using  a  time-integration  technique.  In  this  simulation,  the  first 
thing  computed  is  the  equilibrium  distribution  f0  (x,  k).  Next,  the  voltage  drop  V  is  fixed 
to  0.008  volts,  and  the  electron  distribution/ (x,  k,  t)  is  tracked  for  2000  femtoseconds, 
with  the  initial  condition  set  to/ (x,  k,  0)  =f0  (x,  k).  As  the  distribution  is  being  tracked, 
the  current  density  j(x,  t )  is  computed  from  Eq.  (12).  After  the  2000  femtosecond 
simulation,  the  voltage  drop  V  is  increased  by  0.008  volts,  the  final  electron  distribution 
from  the  previous  time  simulation  is  set  as  the  initial  condition,  and  the  electron 
distribution  is  again  tracked  for  2000  femtoseconds.  This  process  is  repeated  until  a 
predetermined  range  of  voltages  have  been  processed. 
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nx=86,  nk=72  Grid:  Time  Integrator  I -V  Curve 


Fig.  3.  Time-dependent  current- voltage  plot 


The  advantage  of  the  time-dependent  simulation  is  that  if  current  oscillation  develops, 
it  will  be  immediately  noticed  since  the  current  density  is  computed  as  the  electron 
distribution  evolves  in  time.  The  main  disadvantage  of  this  simulation  is  that  it  is 
computationally  intensive.  The  initial  time-integration  technique  took  several  days  to 
compute  the  time-evolution  of  the  electron  distribution  for  the  various  desired  voltage 
drops  for  coarse  grid  studies  (i.e.,  Nx  =  86,  Nt  =  72  or  roughly  6,000  unknowns).  We 
have  improved  upon  the  time-integration  technique  (which  was  previously  a  semi- 
implicit  Euler’s  method,  but  currently  is  an  implicit  Adams  method8)  and  reported  on  this 
in  previous  work9'  I0,  but  we  were  still  restricted  to  coarse  grid  simulation.  Figure  3 
shows  the  current-voltage  plot  generated  by  the  time-dependent  simulation.  This 
represents  the  current  at  the  right-hand  side  (x  =  L)  after  the  electron  distribution  as 
reached  steady-state.  If  no  steady-state  was  reached  but  current  oscillation  occurred,  the 
average  value  of  the  current  was  plotted  instead. 


Current  Density  (A/cm  ) 
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nx=86,  nk=72  Grid 


Fig.  4.  Voltage  drops  that  create  current  oscillation  within  the  RTD 


Figure  4  shows  the  voltage  drops  that  lead  to  current  oscillation  for  the  coarse  grid. 
Here,  V  =  0.248  volts  and  V  =  0.256  volts  are  the  two  voltage  drops  where  current 
oscillation  is  present.  If  one  looks  at  the  current  oscillation  for  V  =  0.248  volts,  the 
period  of  oscillation  appears  to  be  about  400  femtoseconds.  This  corresponds  to  a 
frequency  of  2.5  THz. 

5.  Steady-State  Simulation 

The  second  method  of  simulation  is  the  steady-state  simulation.  We  now  look  directly 
for  steady-state  electron  distributions  as  the  voltage  drop  is  changed.  To  do  this,  we 
solve  a  nonlinear  equation,  given  by  Eq.  (2).  The  advantage  of  this  method  is  that  it  is 
less  computationally  intensive  to  calculate  the  current-voltage  curve,  with  the  steady-state 
simulation  taking  less  than  an  hour  to  compute  the  current-voltage  curve.  If  no  current 
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oscillation  is  occurring  within  the  RTD  for  a  given  voltage  drop  V,  the  time-dependent 
method  will  converge  to  the  steady-state  electron  distribution  which  is  much  more  costly 
than  just  directly  solving  for  the  steady-state  electron  distribution.  Another  advantage  of 
this  method  is  that  fine  grids  can  be  computed.  With  this  method  of  simulation,  we  have 
computed  on  the  Nx  =  1024,  Nk  =  2048  grid,  which  is  roughly  2,000,000  unknowns. 
The  disadvantage  of  this  method  is  that  to  detect  current  oscillation  from  steady-state 
information,  an  eigensolve  must  be  performed  which  is  computationally  expensive.  This 
is  explained  in  the  stability  section. 

To  solve  the  nonlinear  equation  for  the  steady-state  electron  voltages,  we  use 
continuation  methods.  Continuation  methods  solve  nonlinear  equations  which  depend  on 
a  parameter.  For  example,  if  we  are  solving  the  nonlinear  equation  G  (z,  a)  =  0  for  z, 
where  z  is  in  R  m  is  our  state  variable,  and  a,  a  real  number,  is  our  parameter,  continuation 
methods  trace  out  the  solution  curve  z(«),  parameterized  by  a,  so  that  G  (z(«),  a)  =  0.  For 
our  particular  application,  we  are  solving  W (f  V)  =  0,  where/is  the  steady-state  electron 
distribution,  and  the  voltage  drop  V  is  the  parameter.  Numerically,  continuation  methods 
generate  two  sequences:  a  sequence  of  parameters  {V}  (in  our  case,  voltage  drops)  with 
a  corresponding  sequence  of  solutions  {/  '}  (in  our  case,  steady-state  electron  distribution) 
such  that  W  (f1,  V)  =  0.  To  implement  the  continuation  methods  numerically,  we  use  the 
Library  of  Continuation  Algorithms  (LOCA)11. 

LOCA  is  a  part  of  the  Trilinos12  framework,  Sandia  National  Laboratories’  collection 
of  parallel  solver  packages.  Our  implementation  of  LOCA  with  our  RTD  simulator 
makes  use  of  several  other  packages  in  Trilinos,  including  NOX  -  the  nonlinear  solver 
package,  AztecOO  -  the  preconditioned  Krylov  linear  solver  package,  Anasazi  -  the 
eigensolver  package,  and  Epetra  -  a  parallel  data  structure  package.  We  used  a  Newton- 
GMRES  algorithm  to  solve  the  nonlinear  equations,  which  grouped  together  the  NOX 
and  AztecOO  packages,  and  the  Anasazi  eigensolver  to  determine  when  oscillation 
onsets.  To  handle  the  fine  grid  simulations,  we  parallelized  the  RTD  simulator,  which 
involved  using  the  Epetra  package. 

6.  Stability  of  Steady-State  Solutions 

For  the  steady-state,  since  we  are  not  tracking  the  evolution  of  the  current  density  over 
time,  we  need  another  method  to  determine  if  current  oscillation  is  occurring.  Suppose 
we  have  the  nonlinear  ODE  dzldt  =  g  (z)  for  the  state  z  in  R  111  and  a  steady-state  solution 
z  ,  and  we  want  to  determine  the  stability  of  this  steady  state  solution.  One  well-known 
method12  would  be  to  compute  the  eigenvalues  of  the  Jacobian  of  g  (z)  evaluated  at  the 
steady-state  solution  z  ,  denoted  by  g\z  ).  If  all  the  eigenvalues  have  negative  real  part, 
then  the  steady-state  z  is  stable.  If  any  of  the  eigenvalues  have  positive  real  part,  then 
the  steady-state  z  is  unstable.  LOCA  interfaces  to  the  Trilinos  eigensolver,  Anasazi,  to 
compute  the  leading  eigenvalues  of  this  Jacobian  matrix  as  the  continuation  algorithm 
finds  the  steady-state  solutions. 

For  our  application,  we  desire  an  instability  to  develop.  As  the  voltage  drop  V 
changes,  the  eigenvalues  of  the  Jacobian  matrix  at  the  steady-state  electron  distribution 
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will  also  change,  possibly  leading  to  a  change  in  the  stability  of  the  steady-state 
distribution.  This  change  in  stability  is  known  as  a  bifurcation13.  In  particular,  we  want 
to  go  from  a  stable  steady-state  to  oscillatory  behavior  as  the  voltage  drop  V  is  varied. 
This  change  in  stability  is  known  as  a  Hopf  bifurcation.  This  bifurcation  is  characterized 
by  a  complex-conjugate  pair  of  eigenvalues  of  the  Jacobian  whose  real  part  changes  from 
negative  to  positive. 

The  computational  challenge  with  using  the  eigenvalue  method  is  in  detecting  the 
relevant  eigenvalues.  Anasazi  implements  an  iterative  block  Arnoldi  method  for  the 
eigensolve,  which  quickly  locates  the  eigenvalues  with  largest  magnitude.  These 
eigenvalues  are  not  necessarily  the  ones  we  want  to  find,  and  therefore,  the  method  may 
take  a  while  to  locate  the  eigenvalues  of  interest.  One  way  to  handle  this  problem  is  to 
use  a  spectral  transformation  to  move  the  eigenvalues  of  interest  so  that  they  have  larger 
magnitudes.  Let  J  denote  the  Jacobian  matrix.  The  eigenproblem  we  originally  solve  is 
Jz  =  %z.  The  Cayley  transformation  of  the  Jacobian  matrix  is  C  =  (J  -r>I)  '(.J  -  y/),  where 
a  and  y  are  real  numbers.  The  new  eigenproblem  to  solve  is  Cz  =  fiz.  The  eigenvectors 
of  J  and  C  are  the  same,  but  if  o  and  y  are  chosen  so  that  o  >  0  and  y  =  -o,  then  the 
Cayley  transformation  has  the  property14  that  the  Jacobian’s  eigenvalues  which  have 
negative  real  parts  are  mapped  to  the  interior  of  the  unit-circle,  and  the  Jacobian’s 
eigenvalues  which  have  positive  real  parts  are  mapped  to  the  exterior  of  the  unit-circle. 
Since  Anasazi  locates  eigenvalues  of  largest  magnitude  first,  when  solving  the  Cayley 
transformed  eigenproblem,  it  will  quickly  locate  eigenvalues  corresponding  to  the 
unstable  eigenvalues  of  the  Jacobian. 

7.  Steady-state  Numerical  Results 

Figure  5  shows  the  current-voltage  plot  generated  by  the  steady-state  simulation.  This 
represents  the  current  at  the  right-hand  side  (x  =  L)  after  the  electron  distribution  as 
reached  steady-state.  Note  that  these  results  agree  with  what  is  produced  by  the  time- 
dependent  simulation,  except  for  an  additional  branch  connecting  the  higher  part  of  the 
current-voltage  curve  to  the  lower  part  of  the  current  voltage-curve.  This  branch 
represents  solutions  to  the  Wigner-Equations  that  are  unstable  steady-state  electron 
distributions  that  will  not  be  reached  under  time-dependent  simulation. 
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nx=86,  nk=72  Grid:  LOCA  I -V  Curve 


Fig.  5.  Coarse  grid  current- voltage  plot  from  the  steady-state  simulation 


Figure  6  plots  the  complex-conjugate  eigenvalues  of  the  Jacobian  matrix  responsible 
for  the  Hopf  bifurcation  and  onset  of  current  oscillation  on  the  coarse  grid.  The 
eigenvalues  are  tracked,  starting  at  voltage  drop  V  =  0.240  volts,  increasing  the  voltage 
drop  by  0.002  volts,  and  ending  at  V  =  0.256.  At  V  =  0.240  volts,  the  real  parts  of  the 
eigenvalues  are  negative,  so  we  expect  no  oscillation.  This  agrees  with  what  the  time- 
dependent  simulation  produces.  Between  V  =  0.240  volts  and  V  =  0.242  volts,  the 
leading  eigenvalues’  real  parts  become  positive,  indicating  the  onset  of  oscillation.  These 
real  parts  stay  positive  as  V  increases  up  to  0.256  volts,  where  they  are  about  to  turn  back 
negative.  Therefore,  between  V  =  0.242  volts  and  0.256  volts,  we  expect  oscillation,  and 
as  the  voltage  is  further  increased  past  V  =0.256  volts,  we  expect  the  oscillation  to  stop, 
and  a  stable  steady-state  to  return.  All  of  this  agrees  with  what  is  predicted  from  the  time 
dependent  simulation.  Further,  at  V  =  0.248  volts,  the  absolute  value  of  the  imaginary 
parts  of  the  eigenvalues  is  to  =  0.017.  This  corresponds  to  a  period  of  oscillation  of  2ti/co 
=  370  femtoseconds,  which  agrees  closely  with  the  time-dependent  simulation. 


Imaginary  Part 
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Fig.  6.  Complex-conjugate  eigenvalues  that  create  oscillation  on  coarse  grid 


Figure  7  plots  the  coarse  grid  and  fine  grid  current-voltage  plots  generated  by  the 
steady-state  method.  The  fine  grid  current-voltage  curve  is  lower  than  the  coarse  grid 
current-voltage  curve,  and  the  fine  grid  current-voltage  curve  has  more  resonant  features 
than  the  coarse  grid  current-voltage  curve  as  evident  by  the  multiple  peaks  in  the  fine  grid 
result.  These  multiple  peaks  do  not  match  with  experimental  results,  but  we  note  that  for 
the  coarse  grid  result,  the  correlation  length  was  altered  to  match  the  simulated  results  as 
much  as  possible  to  the  experimental  results.  The  correlation  length  will  need  to  be 
readjusted  for  the  fine  grid  simulation  as  part  of  our  ongoing  research. 

8.  Conclusions 

We  have  simulated  electron  transport  within  a  prototypical  nanoscale  semiconductor 
device,  a  resonant  tunneling  diode,  using  the  Wigner-Poisson  equations  in  both  time- 
dependent  and  steady-state  manners.  While  the  time-dependent  method  easily  detects 
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Fig.  7.  Comparison  of  coarse  grid  and  fine  grid  current-voltage  plots 


current  oscillation,  its  computational  cost  makes  this  method  feasible  for  only  coarse 
grids.  In  contrast,  the  steady-state  method  can  handle  fine  grids  and  more  quickly 
generates  the  current-voltage  curve  of  the  RTD  than  the  time-dependent  methods,  but 
detection  of  current  oscillation  requires  a  computationally  intensive  eigensolve.  The 
Cayley  transformation  eases  the  computational  burden  of  the  eigensolve,  and  current 
research  is  focusing  on  finding  a  correlation  length  that  makes  the  current-voltage  curve 
produced  on  the  fine  grids  match  with  experimental  results  and  using  the  Cayley 
transformation  to  locate  the  onset  of  current  oscillation  on  the  fine  grids. 
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